Integration of multimodal cues does not alter mean but reduces variance in avian responses to predators: a systematic review and meta-analysis

1 Setting up

Code

#install.packages("pacman")
#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)

#devtools::install_github("daniel1noble/orchaRd", force = TRUE)

library(tidyverse)
library(here)
library(lme4)
library(orchaRd)
#library(gptstudio)
library(metafor)
library(patchwork)
library(alluvial)
library(ggalluvial)
library(easyalluvial)
library(ape)
library(clubSandwich)
library(emmeans)
library(MuMIn)
library(kableExtra)
library(pander)
# making metafor talk to MuMIn
eval(metafor:::.MuMIn)
# install.packages("pak")
#pak::pak("MichelNivard/gptstudio")

2 Getting data loaded

Code

#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))
#dat_full <- read.csv(here("data/dat_28_06_2023.csv"))
dat_full <- read.csv(here("data/dat_19_07_2023.csv"))

# add phylogenetic tree - only topologies
# TODO? - we could get better tree from birdtree.org
# we can do 50 different trees as in 
# https://academic.oup.com/sysbio/article/68/4/632/5267840
tree_top <- read.tree(here("R/birds_MA.tre"))

# tree with branch lengths
tree <- compute.brlen(tree_top)
#plot(tree)
# turning it into a correlation matrix
cor_tree <- vcv(tree,corr=T)

3 Custom functions

Code

# custom functions

#' Title: Contrast name generator
#'
#' @param name: a vector of character strings
cont_gen <- function(name) {
  combination <- combn(name, 2)
  name_dat <- t(combination)
  names <- paste(name_dat[, 1], name_dat[, 2], sep = "-")
  return(names)
}

#' @title get_pred1: intercept-less model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred1 <- function (model, mod = " ") {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, "")))
  len <- length(name)
  
   if (len != 1) {
        newdata <- matrix(NA, ncol = len, nrow = len)
        for (i in 1:len) {
            pos <- which(model$X[, i] == 1)[[1]]
            newdata[, i] <- model$X[pos, ]
        }
        pred <- metafor::predict.rma(model, newmods = newdata)
    }
    else {
        pred <- metafor::predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title get_pred2: normal model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred2 <- function (model, mod = " ") {
  name <- as.factor(str_replace(row.names(model$beta), 
                                paste0("relevel", "\\(", mod,", ref = name", "\\)"),
                                ""))
  len <- length(name)
  
  if(len != 1){
  newdata <- diag(len)
  pred <- predict.rma(model, intercept = FALSE, newmods = newdata[,-1])
  }
  else {
    pred <- predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title mr_results
#' @description Function to put results of meta-regression and its contrasts
#' @param res1: data frame 1
#' @param res1: data frame 2
mr_results <- function(res1, res2) {
  restuls <-tibble(
    `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)),
    Estimate = c(res1$estimate, res2$estimate),
    `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL),
    `Upper CI  [0.975]` = c(res1$upperCL, res2$upperCL),
    `P value` = c(res1$pval, res2$pval),
    `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR),
    `Upper PI  [0.975]` = c(res1$upperPR, res2$upperPR),
  )
}


#' @title all_models
#' @description Function to take all possible models and get their results
#' @param model: intercept-less model
#' @param mod: the name of a moderator 

all_models <- function(model, mod = " ", type = "homo") {
  
  # getting the level names out
  level_names <- levels(factor(model$data[[mod]]))
  dat2 <- model$data

  VCV1 <- vcalc(vi = dat2$Vd,
             cluster = dat2$SubjectID,
             rho = 0.5)
  #model$data[[mod]] <- factor(model$data[[mod]], ordered = FALSE)
  # meta-regression: contrasts 
  # helper function to run metafor meta-regression
  run_rma1 <- function(name) {
    rma.mv(yi = SMD, 
         V = VCV1, 
         mods = ~ relevel(dat2[[mod]], ref = name), 
         random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat2)
   }

    run_rma2 <- function(name) {
    rma.mv(yi = SMD, 
         V = VCV1, 
         mods = ~ relevel(dat2[[mod]], ref = name), 
         random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~gsub("\"", "", mod) | Obs_ID),
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat2)
   }

  # results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])

  if (type == "homo"){
  model_all <- map(level_names[-length(level_names)], run_rma1)
  } 
  else {
  model_all <- map(level_names[-length(level_names)], run_rma2)
  }
  
  # getting estimates from intercept-less models (means for all the groups)
  res1 <- get_pred1(model, mod = mod)
  
  # getting estiamtes from all contrast models
  res2_pre <- map(model_all, ~ get_pred2(.x, mod = mod))
  
  # a list of the numbers to take out unnecessary contrasts
  contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1))
  res2 <- map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) 
  # creating a table
  res_tab <- mr_results(res1, res2) %>% 
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")
  
  # results
  res_tab

}

# test 
#all_models(mod1, mod = "Treat_mod")
#all_models(mod1s, mod = "Treat_mod")

4 Preparing data set

Code

# function for calculating variance
Vd_func <- function(d, n1, n2, design, r = 0.5){
  # independent design
  if(design == "among"){
    var <- (n1 + n2) / (n1*n2) + d^2 / (2 * (n1 + n2 - 2)) # variance
  } else { # dependent design
    var <- 2*(1-r) / n1 + d^2 / (2*(n1 - 1)) # variance
  }
  var # return variance
}

# getting Hedges' g - get small size bias corrected effect size
dat_full$SMD <- dat_full$d / (1 - 3/(4 * (dat_full$NTreat + dat_full$Ncontrol) - 9))

# flipping d 
dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection


# calucating Vd
dat_full$Vd <- with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))


# extra useful function
# function for getting mean and sd from median, quartiles and sample size
# get_mean_sd <- function(median, q1, q3, n){
#   sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd
#   mean <- (median + q1 + q3)/3 # mean
#   c(mean, sd)
# }


# observation id
dat_full$Obs_ID <- 1:nrow(dat_full)
dat_full$Phylo <- gsub(" ", "_", dat_full$FocalSpL)

# filtering very large variance and also very small sample size
dat_int <- dat_full %>% filter(Vd < 10 & Ncontrol > 2 & NTreat > 2)

#dim(dat_full)
#dim(dat_int)


# sorting out modality stuff
# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) 

dat_int %>% mutate(Treat_mod = case_when(Treatment == "A" ~ "A",
                                          Treatment == "AV" ~ "AV",
                                          Treatment == "AVG" ~ "AV",
                                          Treatment == "AVM" ~ "AV",
                                          Treatment == "L" ~ "AVO",
                                          Treatment == "O" ~ "O",
                                          Treatment == "OV" ~ "OV",
                                          Treatment == "V" ~ "V",
                                          Treatment == "VG" ~ "V",
                                          Treatment == "VM" ~ "V",
                                          Treatment == "VP" ~ "V"),
                    # into how many
                    Treat_No = case_when(Treatment == "A" ~ 1,
                                         Treatment == "AV" ~ 2,
                                         Treatment == "AVG" ~ 2,
                                         Treatment == "AVM" ~ 2,
                                         Treatment == "L" ~ 3,
                                         Treatment == "O" ~ 1,
                                         Treatment == "OV" ~ 2,
                                         Treatment == "V" ~ 1,
                                         Treatment == "VG" ~ 1,
                                         Treatment == "VM" ~ 1,
                                         Treatment == "VP" ~ 1),
                    # des it have some add-ons
                    Add_on = case_when(Treatment == "A" ~ "No",
                                         Treatment == "AV" ~ "No",
                                         Treatment == "AVG" ~ "Yes",
                                         Treatment == "AVM" ~ "Yes",
                                         Treatment == "L" ~ "No",
                                         Treatment == "O" ~ "No",
                                         Treatment == "OV" ~ "No",
                                         Treatment == "V" ~ "No",
                                         Treatment == "VG" ~ "Yes",
                                         Treatment == "VM" ~ "Yes",
                                         Treatment == "VP" ~ "Yes"),

                      ) -> dat

# creating data just for A, V, and AV 
dat_short <- dat %>% filter(Treat_mod == "A" | Treat_mod == "V" | Treat_mod == "AV")

# for add-on, we only need V and AV
dat_short_add <- dat %>% filter(Treat_mod == "AV" | Treat_mod == "V")


dat <- dat %>%
  mutate_if(is.character, as.factor)

# reordering factors for better visualisation

# dat$Treat_mod <- factor(dat$Treat_mod, levels = c("A", "V", "AV", "O", "OV", "AVO"))

5 Exploratory visualization

For Treat_mod (Treatment), we will only visualise A, V, and AV as O $r $, OV, and AVO are much rarer. But for Type (Trait type), we will use all data.

Treatment vs. Trait types

Code

# reordering dat_shrot$Treat_mod for better visualisation
dat_short$Treat_mod <- factor(dat_short$Treat_mod, levels = c("A", "V", "AV"))

# Treatment vs Type
dat_short %>% group_by(Treat_mod, Type) %>%
  summarise(n = n()) -> tab

#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Type)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and trait type")

Treatment vs. Treatment duration

Code

# tuning Treatment duration into a binary variable
dat_short %>% mutate(TDration = case_when(duration_days < 1 ~ "< 1 day",
                                          duration_days >= 1 ~ "> 1 day")) -> dat_short

# reformatting data
dat_short %>% group_by(Treat_mod, TDration) %>%
  summarise(n = n()) -> tab

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = TDration)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and duration of treatment")

Treatment vs. Sex

Code

 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Sex) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Sex)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs. Predator type

Code

# Treat_mod vs Design
dat_short %>% filter(PredTo != "") %>%
  group_by(Treat_mod, PredTo) %>%
  summarise(n = n()) -> tab

tab$PredTo <- factor(tab$PredTo, levels = c("A", "N", "B"), labels = c("Adult", "Nestling", "Both"))

#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = PredTo)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs. Design

Code

 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Design) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Design)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs. Season

Code

 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Season) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Season)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs. Setting

Code

 # Treat_mod vs Design

dat_short %>% group_by(Treat_mod, Setting) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Setting)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and study setting")

Treatment vs Control type

Code

# Treat_mod vs ControlType

dat_short %>% filter(ControlType != "mix") %>% group_by(Treat_mod, ControlType) %>%
  summarise(n = n()) -> tab
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = ControlType)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 6, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and control type")

Trait type vs. Treatment duration

Code

# Type vs duration_days

# turn duration_days into a binary factor: less than 1 and more than 1

dat %>% mutate(TDration = case_when(duration_days < 1 ~ "< 1 day",
                                          duration_days >= 1 ~ "> 1 day")) -> dat


dat %>% group_by(Type, TDration) %>%
  summarise(n = n()) -> tab2

# using ggaruvial
ggplot(tab2,
       aes(y = n,
           axis1 = Type,
           axis2 = TDration)) +
  geom_alluvium(aes(fill = Type)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Trait type and duration of treatment")

Trait type vs. Sex

Code

# Type vs duration_days

# turn duration_days into a binary factor: less than 1 and more than 1

dat %>% mutate(TDration = case_when(duration_days < 1 ~ "< 1 day",
                                          duration_days >= 1 ~ "> 1 day")) -> dat

dat %>% group_by(Type, Sex) %>%
  summarise(n = n()) -> tab

# reordering
tab$Sex <- factor(tab$Sex, levels = c("F", "M", "both"), labels = c("Female", "Male", "Both"))


# using ggaruvial
ggplot(tab,
       aes(y = n,
           axis1 = Type,
           axis2 = Sex)) +
  geom_alluvium(aes(fill = Type)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 4, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Trait type and duration of treatment")

6 Meta-analysis

All random effects

Code

##| warning: false
# VCV matrix

VCV <- vcalc(vi = dat$Vd,
             cluster = dat$SubjectID,
             rho = 0.5)

mod0 <- rma.mv(yi = SMD,
       V = VCV, 
       random = list(#~1 | Phylo,
                     ~1 | FocalSpL,
                     ~1 | RecNo,
                     ~1 | SubjectID, # incoprated as VCV
                     ~1 | Obs_ID),
      # R = list(Phylo = cor_tree),
       test = "t",
       method = "REML", 
       sparse = TRUE,
       data = dat)

summary(mod0)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-894.6453  1789.2906  1799.2906  1821.5901  1799.3853   

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.0111  0.1053     90     no   FocalSpL 
sigma^2.2  0.1450  0.3808    116     no      RecNo 
sigma^2.3  0.0000  0.0001    153     no  SubjectID 
sigma^2.4  0.6872  0.8290    640     no     Obs_ID 

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.4139  0.0653  6.3401  639  <.0001  0.2857  0.5421  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# TODO - think about whether we add this or not
robust(mod0, cluster = dat$SubjectID)

Multivariate Meta-Analysis Model (k = 640; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.0111  0.1053     90     no   FocalSpL 
sigma^2.2  0.1450  0.3808    116     no      RecNo 
sigma^2.3  0.0000  0.0001    153     no  SubjectID 
sigma^2.4  0.6872  0.8290    640     no     Obs_ID 

Number of estimates:   640
Number of clusters:    153
Estimates per cluster: 1-29 (mean: 4.18, median: 2)

Model Results:

estimate      se¹    tval¹   df¹    pval¹   ci.lb¹   ci.ub¹      
  0.4139  0.0540   7.6691   152   <.0001   0.3073   0.5206   *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1) results based on cluster-robust inference (var-cov estimator: CR1,
   approx t-test and confidence interval, df: residual method)

Code

round(i2_ml(mod0), 2)

    I2_Total  I2_FocalSpL     I2_RecNo I2_SubjectID    I2_Obs_ID 
       92.81         1.22        15.96         0.00        75.63 

Code

orchard_plot(mod0,
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Reduced model

Code

# reduced model

mod0r <- rma.mv(yi = SMD,
       V = VCV, 
       random = list(~1 | FocalSpL,
                     ~1 | RecNo,
                     ~1 | Obs_ID),
       test = "t",
       method = "REML", 
       sparse = TRUE,
       data = dat)

summary(mod0r)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-894.6453  1789.2906  1797.2906  1815.1302  1797.3537   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0111  0.1052     90     no  FocalSpL 
sigma^2.2  0.1450  0.3808    116     no     RecNo 
sigma^2.3  0.6872  0.8290    640     no    Obs_ID 

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.4139  0.0653  6.3401  639  <.0001  0.2857  0.5421  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(i2_ml(mod0r), 2)

   I2_Total I2_FocalSpL    I2_RecNo   I2_Obs_ID 
      92.81        1.22       15.96       75.63 

Code

# comparing two models
#anova(mod0, mod0r)

7 Meta-regression: uni-moderator (mostly)

Treatmeant with all data

Code

## Treatment - A, V, AV etc 
#dat$Phylo <- as.character(dat$Phylo)
#match(dat$Phylo, colnames(cor_tree))
#match(colnames(cor_tree), dat$Phylo)

mod1 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Treat_mod - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod1)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-887.7086  1775.4172  1793.4172  1833.4856  1793.7056   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0074  0.0862     90     no  FocalSpL 
sigma^2.2  0.1563  0.3954    116     no     RecNo 
sigma^2.3  0.6907  0.8311    640     no    Obs_ID 

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 634) = 6.8763, p-val < .0001

Model Results:

              estimate      se    tval   df    pval    ci.lb   ci.ub      
Treat_modA      0.4019  0.0971  4.1381  634  <.0001   0.2112  0.5926  *** 
Treat_modAV     0.4723  0.1345  3.5128  634  0.0005   0.2083  0.7363  *** 
Treat_modAVO    0.5955  0.4487  1.3270  634  0.1850  -0.2857  1.4767      
Treat_modO      0.2577  0.2930  0.8797  634  0.3793  -0.3176  0.8330      
Treat_modOV     0.0185  0.3829  0.0482  634  0.9615  -0.7334  0.7703      
Treat_modV      0.4324  0.1043  4.1442  634  <.0001   0.2275  0.6373  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod1)*100, 2)

   R2_marginal R2_conditional 
          0.64          19.68 

Code

orchard_plot(mod1, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# result table
all_models(mod1, mod = "Treat_mod")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
A	0.402	0.211	0.593	0.000	-1.423	2.227
AV	0.472	0.208	0.736	0.000	-1.362	2.307
AVO	0.595	-0.286	1.477	0.185	-1.422	2.613
O	0.258	-0.318	0.833	0.379	-1.646	2.162
OV	0.018	-0.733	0.770	0.962	-1.946	1.983
V	0.432	0.228	0.637	0.000	-1.394	2.259
A-AV	0.070	-0.248	0.389	0.664	-1.772	1.913
A-AVO	0.194	-0.706	1.093	0.673	-1.832	2.219
A-O	-0.144	-0.748	0.460	0.639	-2.057	1.769
A-OV	-0.383	-1.158	0.392	0.332	-2.357	1.590
A-V	0.031	-0.240	0.301	0.825	-1.805	1.866
AV-AVO	0.123	-0.797	1.043	0.793	-1.912	2.158
AV-O	-0.215	-0.847	0.418	0.506	-2.137	1.708
AV-OV	-0.454	-1.250	0.343	0.264	-2.436	1.528
AV-V	-0.040	-0.340	0.260	0.794	-1.880	1.800
AVO-O	-0.338	-1.390	0.714	0.529	-2.436	1.760
AVO-OV	-0.577	-1.711	0.557	0.318	-2.717	1.563
AVO-V	-0.163	-1.067	0.741	0.723	-2.191	1.865
O-OV	-0.239	-1.186	0.707	0.620	-2.286	1.808
O-V	0.175	-0.434	0.784	0.574	-1.740	2.089
OV-V	0.414	-0.364	1.192	0.296	-1.561	2.389

Treatment with Audio and Visual data

7.1 Homoscedastic model

Code

VCVs <- vcalc(vi = dat_short$Vd,
             cluster = dat_short$SubjectID,
             rho = 0.5)


mod1s <- rma.mv(yi = SMD,
                   V = VCVs,
                   random = list(~1|FocalSpL,
                                 ~1 | RecNo,
                                 ~1 | Obs_ID),
                   mod = ~ Treat_mod - 1, 
                   test = "t",
                   method = "REML", 
                   sparse = TRUE,
                   data = dat_short)

summary(mod1s)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-847.7054  1695.4108  1707.4108  1733.7623  1707.5532   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0087  0.0932     81     no  FocalSpL 
sigma^2.2  0.1700  0.4123    104     no     RecNo 
sigma^2.3  0.7138  0.8449    600     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 597) = 12.6006, p-val < .0001

Model Results:

             estimate      se    tval   df    pval   ci.lb   ci.ub      
Treat_modA     0.4055  0.0996  4.0715  597  <.0001  0.2099  0.6012  *** 
Treat_modV     0.4386  0.1069  4.1041  597  <.0001  0.2287  0.6485  *** 
Treat_modAV    0.4782  0.1379  3.4687  597  0.0006  0.2074  0.7490  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod1s)*100, 2)

   R2_marginal R2_conditional 
          0.08          20.09 

Code

orchard_plot(mod1s, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# result table
all_models(mod1s, mod = "Treat_mod")

Fixed effect	Estimate	Lower CI [0.025]	Upper CI [0.975]	P value	Lower PI [0.025]	Upper PI [0.975]
A	0.406	0.210	0.601	0.000	-1.460	2.271
V	0.439	0.229	0.648	0.000	-1.429	2.306
AV	0.478	0.207	0.749	0.001	-1.397	2.353
A-V	0.033	-0.244	0.310	0.815	-1.843	1.909
A-AV	0.073	-0.254	0.399	0.662	-1.811	1.957
V-AV	0.040	-0.267	0.346	0.800	-1.841	1.920

7.2 Heteroscedastic model

Code

# modeling heteroscedasticity
mod1s2 <- rma.mv(yi = SMD, 
                V = VCVs, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~ Treat_mod | Obs_ID), 
                mod = ~ Treat_mod, 
                test = "t",
                struct = "DIAG",
                method = "REML", 
                sparse = TRUE,
                data = dat_short)

summary(mod1s2)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-838.2594  1676.5188  1692.5188  1727.6541  1692.7637   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0086  0.0927     81     no  FocalSpL 
sigma^2.2  0.1764  0.4200    104     no     RecNo 

outer factor: Obs_ID    (nlvls = 600)
inner factor: Treat_mod (nlvls = 3)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.7479  0.8648    302     no      A 
tau^2.2    0.8583  0.9264    190     no      V 
tau^2.3    0.3460  0.5882    108     no     AV 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 597) = 0.0343, p-val = 0.9662

Model Results:

             estimate      se    tval   df    pval    ci.lb   ci.ub      
intrcpt        0.4100  0.1013  4.0495  597  <.0001   0.2112  0.6089  *** 
Treat_modV     0.0300  0.1452  0.2068  597  0.8363  -0.2552  0.3152      
Treat_modAV    0.0382  0.1550  0.2466  597  0.8053  -0.2662  0.3426      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mod1s2, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

7.3 Comparing models

Code

# comparision models
anova(mod1s, mod1s2)

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1692.5188 1727.6541 1692.7637 -838.2594                NA 
Reduced  6 1707.4108 1733.7623 1707.5532 -847.7054 18.8921 <.0001 NA 

Treatment with some additions

Code

# the effect of additions
# this is a part of sensitivity analysis

VCVs2 <- vcalc(vi = dat_short_add$Vd,
             cluster = dat_short_add$SubjectID,
             rho = 0.5)

mod5 <- rma.mv(yi = SMD, 
                V = VCVs2, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~ Treat_mod | Obs_ID), 
                mod = ~ Treat_mod*Add_on - 1,
                test = "t",
                struct = "DIAG",
                method = "REML", 
                sparse = TRUE,
                data = dat_short_add)

summary(mod5)

Multivariate Meta-Analysis Model (k = 298; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-414.2668   828.5336   844.5336   874.0022   845.0388   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     39     no  FocalSpL 
sigma^2.2  0.1561  0.3950     58     no     RecNo 

outer factor: Obs_ID    (nlvls = 298)
inner factor: Treat_mod (nlvls = 2)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.3485  0.5904    108     no     AV 
tau^2.2    0.8462  0.9199    190     no      V 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 294) = 6.0475, p-val = 0.0001

Model Results:

                      estimate      se     tval   df    pval    ci.lb   ci.ub 
Treat_modAV             0.4064  0.1261   3.2223  294  0.0014   0.1582  0.6546 
Treat_modV              0.4344  0.1166   3.7269  294  0.0002   0.2050  0.6639 
Add_onYes               0.1193  0.2685   0.4442  294  0.6573  -0.4092  0.6477 
Treat_modV:Add_onYes   -0.2433  0.3460  -0.7033  294  0.4824  -0.9242  0.4376 
                          
Treat_modAV            ** 
Treat_modV            *** 
Add_onYes                 
Treat_modV:Add_onYes      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Treatment as an ordinal variable

Code

# testing the number of stimuli

mod4 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~1 | Obs_ID), 
               mod = ~ Treat_No, 
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod4)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.4118  1786.8236  1796.8236  1819.1153  1796.9185   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0117  0.1081     90     no  FocalSpL 
sigma^2.2  0.1469  0.3833    116     no     RecNo 
sigma^2.3  0.6875  0.8292    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 0.0952, p-val = 0.7577

Model Results:

          estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt     0.3696  0.1603  2.3057  638  0.0214   0.0548  0.6843  * 
Treat_No    0.0366  0.1185  0.3086  638  0.7577  -0.1962  0.2694    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(mod4,
             mod = "Treat_No",
             group = "RecNo",
             xlab = "The number of simuli",
             g = TRUE)

Trait type

Code

# Type of responses
mod2 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Type - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod2)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-886.7295  1773.4590  1785.4590  1812.1996  1785.5923   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0002     90     no  FocalSpL 
sigma^2.2  0.1322  0.3635    116     no     RecNo 
sigma^2.3  0.6902  0.8308    640     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 18.0966, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
TypeBehaviour      0.4866  0.0677  7.1917  637  <.0001   0.3537  0.6194  *** 
TypeLifeHistory    0.3089  0.1176  2.6263  637  0.0088   0.0779  0.5398   ** 
TypePhysiology     0.0284  0.1285  0.2207  637  0.8254  -0.2240  0.2807      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod2c <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Type,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod2c)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-886.7295  1773.4590  1785.4590  1812.1996  1785.5923   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0002     90     no  FocalSpL 
sigma^2.2  0.1322  0.3635    116     no     RecNo 
sigma^2.3  0.6902  0.8308    640     no    Obs_ID 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 637) = 5.9194, p-val = 0.0028

Model Results:

                 estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt            0.4866  0.0677   7.1917  637  <.0001   0.3537   0.6194  *** 
TypeLifeHistory   -0.1777  0.1222  -1.4549  637  0.1462  -0.4176   0.0622      
TypePhysiology    -0.4582  0.1350  -3.3950  637  0.0007  -0.7232  -0.1932  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod2d <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ relevel(Type, ref = "LifeHistory"),
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod2d)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-886.7295  1773.4590  1785.4590  1812.1996  1785.5923   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0002     90     no  FocalSpL 
sigma^2.2  0.1322  0.3635    116     no     RecNo 
sigma^2.3  0.6902  0.8308    640     no    Obs_ID 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 637) = 5.9194, p-val = 0.0028

Model Results:

                                              estimate      se     tval   df 
intrcpt                                         0.3089  0.1176   2.6263  637 
relevel(Type, ref = "LifeHistory")Behaviour     0.1777  0.1222   1.4549  637 
relevel(Type, ref = "LifeHistory")Physiology   -0.2805  0.1557  -1.8017  637 
                                                pval    ci.lb   ci.ub     
intrcpt                                       0.0088   0.0779  0.5398  ** 
relevel(Type, ref = "LifeHistory")Behaviour   0.1462  -0.0622  0.4176     
relevel(Type, ref = "LifeHistory")Physiology  0.0721  -0.5862  0.0252   . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mod2,
             mod = "Type",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

round(r2_ml(mod2)*100, 2)

   R2_marginal R2_conditional 
          3.15          18.71 

Code

# heteroscadasticity model
mod2b <- rma.mv(yi = SMD, 
               V = VCV, 
               mod = ~ Type - 1, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~Type | Obs_ID), 
               struct = "DIAG",
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod2b)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-870.7865  1741.5731  1757.5731  1793.2272  1757.8024   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0274  0.1656     90     no  FocalSpL 
sigma^2.2  0.1296  0.3600    116     no     RecNo 

outer factor: Obs_ID (nlvls = 640)
inner factor: Type   (nlvls = 3)

            estim    sqrt  k.lvl  fixed        level 
tau^2.1    0.6577  0.8110    434     no    Behaviour 
tau^2.2    0.3176  0.5636    112     no  LifeHistory 
tau^2.3    1.2391  1.1131     94     no   Physiology 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 17.4552, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
TypeBehaviour      0.5026  0.0720  6.9850  637  <.0001   0.3613  0.6439  *** 
TypeLifeHistory    0.3531  0.1047  3.3712  637  0.0008   0.1474  0.5587  *** 
TypePhysiology     0.0120  0.1549  0.0773  637  0.9384  -0.2922  0.3161      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# make other hetero
mod2e <- rma.mv(yi = SMD,
               V = VCV,
               mod = ~ Type, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~Type | Obs_ID), 
               struct = "DIAG",
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod2e)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-870.7865  1741.5731  1757.5731  1793.2272  1757.8024   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0274  0.1656     90     no  FocalSpL 
sigma^2.2  0.1296  0.3600    116     no     RecNo 

outer factor: Obs_ID (nlvls = 640)
inner factor: Type   (nlvls = 3)

            estim    sqrt  k.lvl  fixed        level 
tau^2.1    0.6577  0.8110    434     no    Behaviour 
tau^2.2    0.3176  0.5636    112     no  LifeHistory 
tau^2.3    1.2391  1.1131     94     no   Physiology 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 637) = 5.0637, p-val = 0.0066

Model Results:

                 estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt            0.5026  0.0720   6.9850  637  <.0001   0.3613   0.6439  *** 
TypeLifeHistory   -0.1495  0.1037  -1.4414  637  0.1500  -0.3533   0.0542      
TypePhysiology    -0.4906  0.1567  -3.1304  637  0.0018  -0.7984  -0.1829   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod2f <- rma.mv(yi = SMD,
               V = VCV,
               mod = ~ relevel(Type, ref = "LifeHistory"), 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~Type | Obs_ID), 
               struct = "DIAG",
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod2f)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-870.7865  1741.5731  1757.5731  1793.2272  1757.8024   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0274  0.1656     90     no  FocalSpL 
sigma^2.2  0.1296  0.3600    116     no     RecNo 

outer factor: Obs_ID (nlvls = 640)
inner factor: Type   (nlvls = 3)

            estim    sqrt  k.lvl  fixed        level 
tau^2.1    0.6577  0.8110    434     no    Behaviour 
tau^2.2    0.3176  0.5636    112     no  LifeHistory 
tau^2.3    1.2391  1.1131     94     no   Physiology 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 637) = 5.0637, p-val = 0.0066

Model Results:

                                              estimate      se     tval   df 
intrcpt                                         0.3531  0.1047   3.3712  637 
relevel(Type, ref = "LifeHistory")Behaviour     0.1495  0.1037   1.4414  637 
relevel(Type, ref = "LifeHistory")Physiology   -0.3411  0.1614  -2.1131  637 
                                                pval    ci.lb    ci.ub      
intrcpt                                       0.0008   0.1474   0.5587  *** 
relevel(Type, ref = "LifeHistory")Behaviour   0.1500  -0.0542   0.3533      
relevel(Type, ref = "LifeHistory")Physiology  0.0350  -0.6581  -0.0241    * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# heteroscadasticity model better than the homoscedasticity model
# note LifeHistory has lowest variation but this may be expected? 
# as it is less flexiable (e.g. the number of eggs?)
anova(mod2, mod2b)

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1757.5731 1793.2272 1757.8024 -870.7865                NA 
Reduced  6 1785.4590 1812.1996 1785.5923 -886.7295 31.8859 <.0001 NA 

Code

orchard_plot(mod2b, 
             mod = "Type",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

7.4 Trait categories

Code

# Category of responses

mod3 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~1 | Obs_ID), 
               mod = ~ Category - 1, 
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod3)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-879.5908  1759.1816  1781.1816  1830.1194  1781.6074   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     90     no  FocalSpL 
sigma^2.2  0.1389  0.3727    116     no     RecNo 
sigma^2.3  0.6912  0.8314    640     no    Obs_ID 

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 632) = 7.2433, p-val < .0001

Model Results:

                          estimate      se    tval   df    pval    ci.lb 
CategoryAntiPredator        0.4772  0.0924  5.1669  632  <.0001   0.2959 
CategoryCondition           0.0027  0.1394  0.0191  632  0.9848  -0.2710 
CategoryCostlyBehaviours    0.4567  0.1784  2.5596  632  0.0107   0.1063 
CategoryHormones            0.1717  0.2914  0.5891  632  0.5560  -0.4006 
CategoryIntake              0.7554  0.1623  4.6539  632  <.0001   0.4366 
CategoryParentalCare        0.4103  0.1076  3.8133  632  0.0002   0.1990 
CategoryPhenology           0.1927  0.1799  1.0711  632  0.2845  -0.1606 
Categoryreproduction        0.3366  0.1287  2.6159  632  0.0091   0.0839 
                           ci.ub      
CategoryAntiPredator      0.6586  *** 
CategoryCondition         0.2763      
CategoryCostlyBehaviours  0.8071    * 
CategoryHormones          0.7439      
CategoryIntake            1.0741  *** 
CategoryParentalCare      0.6216  *** 
CategoryPhenology         0.5459      
Categoryreproduction      0.5892   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod3)*100, 2)

   R2_marginal R2_conditional 
          4.12          20.17 

Code

orchard_plot(mod3, 
             mod = "Category",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             angle = 45,
             branch.size = 3)

7.5 Predactor guild

Code

# Predactor guild
# quite heterogeneous
# TODO this could be in random effects - think abou thtis a bit later
mod6 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ PredGuild - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod6)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-886.0833  1772.1665  1790.1665  1830.2350  1790.4550   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0231  0.1519     90     no  FocalSpL 
sigma^2.2  0.1418  0.3765    116     no     RecNo 
sigma^2.3  0.6884  0.8297    640     no    Obs_ID 

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 634) = 7.4143, p-val < .0001

Model Results:

                  estimate      se     tval   df    pval    ci.lb   ci.ub      
PredGuildBird       0.4848  0.0790   6.1335  634  <.0001   0.3296  0.6400  *** 
PredGuildFish      -0.1871  0.9107  -0.2054  634  0.8373  -1.9755  1.6013      
PredGuildMammal     0.2988  0.1305   2.2904  634  0.0223   0.0426  0.5551    * 
PredGuildmixed      0.1133  0.2515   0.4505  634  0.6525  -0.3805  0.6072      
PredGuildNS         0.6840  0.3025   2.2615  634  0.0241   0.0901  1.2779    * 
PredGuildReptile    0.3365  0.2838   1.1857  634  0.2362  -0.2208  0.8939      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod6)*100, 2)

   R2_marginal R2_conditional 
          2.20          21.09 

Code

orchard_plot(mod6, 
             mod = "PredGuild",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Setting

Code

# Setting

mod7 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Setting - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod7)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-891.0585  1782.1169  1794.1169  1820.8575  1794.2502   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0172  0.1311     90     no  FocalSpL 
sigma^2.2  0.1378  0.3712    116     no     RecNo 
sigma^2.3  0.6868  0.8288    640     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 14.2658, p-val < .0001

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub 
SettingField           0.4245  0.0721   5.8846  637  <.0001   0.2829  0.5662 
SettingLab             0.4764  0.1575   3.0247  637  0.0026   0.1671  0.7857 
SettingSemi-natural   -0.2054  0.4204  -0.4886  637  0.6253  -1.0308  0.6201 
                         
SettingField         *** 
SettingLab            ** 
SettingSemi-natural      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod7)*100, 2)

   R2_marginal R2_conditional 
          0.70          18.98 

Code

orchard_plot(mod7, 
             mod = "Setting",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Season

Code

# Season

mod8 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Season - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

mod8b <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Season,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod8b)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-892.1745  1784.3489  1794.3489  1816.6406  1794.4439   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0187  0.1368     90     no  FocalSpL 
sigma^2.2  0.1390  0.3728    116     no     RecNo 
sigma^2.3  0.6875  0.8291    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 2.3205, p-val = 0.1282

Model Results:

                    estimate      se    tval   df    pval    ci.lb   ci.ub      
intrcpt               0.3589  0.0783  4.5831  638  <.0001   0.2051  0.5127  *** 
Seasonnon-breeding    0.2007  0.1318  1.5233  638  0.1282  -0.0580  0.4595      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod8)*100, 2)

   R2_marginal R2_conditional 
          0.95          19.44 

Code

orchard_plot(mod8,
             mod = "Season",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Design

Code

# Design
mod9 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Design - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod9)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.3891  1786.7783  1796.7783  1819.0700  1796.8732   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0103  0.1016     90     no  FocalSpL 
sigma^2.2  0.1481  0.3848    116     no     RecNo 
sigma^2.3  0.6879  0.8294    640     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 638) = 19.9969, p-val < .0001

Model Results:

              estimate      se    tval   df    pval   ci.lb   ci.ub      
Designamong     0.4061  0.0874  4.6483  638  <.0001  0.2345  0.5777  *** 
Designwithin    0.4197  0.0839  5.0006  638  <.0001  0.2549  0.5844  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod9b <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Design,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod9b)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.3891  1786.7783  1796.7783  1819.0700  1796.8732   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0103  0.1016     90     no  FocalSpL 
sigma^2.2  0.1481  0.3848    116     no     RecNo 
sigma^2.3  0.6879  0.8294    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 0.0150, p-val = 0.9026

Model Results:

              estimate      se    tval   df    pval    ci.lb   ci.ub      
intrcpt         0.4061  0.0874  4.6483  638  <.0001   0.2345  0.5777  *** 
Designwithin    0.0135  0.1106  0.1225  638  0.9026  -0.2037  0.2308      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod9)*100, 2)

   R2_marginal R2_conditional 
          0.01          18.72 

Code

orchard_plot(mod9,
             mod = "Design",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Response period

Code

# Response period
mod10 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ResponsePeriod - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod10)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-892.5749  1785.1498  1795.1498  1817.4415  1795.2448   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0138  0.1176     90     no  FocalSpL 
sigma^2.2  0.1435  0.3788    116     no     RecNo 
sigma^2.3  0.6883  0.8296    640     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 638) = 20.8414, p-val < .0001

Model Results:

                      estimate      se    tval   df    pval   ci.lb   ci.ub 
ResponsePeriodafter     0.3500  0.0851  4.1108  638  <.0001  0.1828  0.5171 
ResponsePeriodduring    0.4938  0.0903  5.4708  638  <.0001  0.3166  0.6711 
                          
ResponsePeriodafter   *** 
ResponsePeriodduring  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod10)*100, 2)

   R2_marginal R2_conditional 
          0.61          19.10 

Code

orchard_plot(mod10,
             mod = "ResponsePeriod",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Control type

Code

# control type
mod11 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ControlType - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod11)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-891.3557  1782.7115  1796.7115  1827.8979  1796.8898   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0080  0.0893     90     no  FocalSpL 
sigma^2.2  0.1504  0.3879    116     no     RecNo 
sigma^2.3  0.6901  0.8307    640     no    Obs_ID 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 636) = 10.0508, p-val < .0001

Model Results:

                        estimate      se    tval   df    pval    ci.lb   ci.ub 
ControlTypeBlank          0.4424  0.0947  4.6709  636  <.0001   0.2564  0.6283 
ControlTypedisturbance    0.4170  0.1452  2.8710  636  0.0042   0.1318  0.7022 
ControlTypemix            0.4533  0.9429  0.4807  636  0.6309  -1.3983  2.3049 
ControlTypeNonPred        0.3895  0.0810  4.8106  636  <.0001   0.2305  0.5485 
                            
ControlTypeBlank        *** 
ControlTypedisturbance   ** 
ControlTypemix              
ControlTypeNonPred      *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod11)*100, 2)

   R2_marginal R2_conditional 
          0.07          18.73 

Code

orchard_plot(mod11,
             mod = "ControlType",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Sex

Code

# sex
# TODO - this could be interesting
# what is in males and females
mod12 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Sex - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod12)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-890.5710  1781.1420  1793.1420  1819.8826  1793.2753   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0099  0.0993     90     no  FocalSpL 
sigma^2.2  0.1354  0.3680    116     no     RecNo 
sigma^2.3  0.6906  0.8310    640     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 637) = 15.0441, p-val < .0001

Model Results:

         estimate      se    tval   df    pval   ci.lb   ci.ub      
Sexboth    0.4706  0.0758  6.2110  637  <.0001  0.3218  0.6193  *** 
SexF       0.2579  0.1011  2.5501  637  0.0110  0.0593  0.4565    * 
SexM       0.4309  0.1325  3.2515  637  0.0012  0.1707  0.6912   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod12)*100, 2)

   R2_marginal R2_conditional 
          0.95          18.17 

Code

orchard_plot(mod12,
             mod = "Sex",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# shoter data for just males and females
# hetero but no sex effect here
dat_sex <- dat %>% filter(Sex != "both")

VCV3 <- vcalc(vi = dat_sex$Vd,
             cluster = dat_sex$SubjectID,
             rho = 0.5)

mod12a <- rma.mv(yi = SMD, 
                 V = VCV3, 
                 mod = ~ Sex, 
                 random = list(~1|FocalSpL , 
                               ~1 | RecNo, 
                               ~1 | Obs_ID), 
                 #struct = "DIAG",
                 test = "t",
                 method = "REML", 
                 sparse = TRUE,
                 data = dat_sex)

mod12b <- rma.mv(yi = SMD, 
                V = VCV3, 
                mod = ~ Sex, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~Sex | Obs_ID), 
                struct = "DIAG",
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat_sex)

summary(mod12b)

Multivariate Meta-Analysis Model (k = 242; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-260.7712   521.5424   533.5424   554.4262   533.9029   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0983  0.3136     36     no  FocalSpL 
sigma^2.2  0.0000  0.0002     49     no     RecNo 

outer factor: Obs_ID (nlvls = 242)
inner factor: Sex    (nlvls = 2)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.3015  0.5491    158     no      F 
tau^2.2    0.6644  0.8151     84     no      M 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 240) = 2.7884, p-val = 0.0963

Model Results:

         estimate      se     tval   df    pval    ci.lb   ci.ub      
intrcpt    0.4995  0.1287   3.8826  240  0.0001   0.2461  0.7529  *** 
SexF      -0.2049  0.1227  -1.6698  240  0.0963  -0.4466  0.0368    . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

anova(mod12a, mod12b)

        df      AIC      BIC     AICc    logLik     LRT   pval QE 
Full     6 533.5424 554.4262 533.9029 -260.7712                NA 
Reduced  5 542.3039 559.7071 542.5604 -266.1520 10.7616 0.0010 NA 

Code

orchard_plot(mod12b,
             mod = "Sex",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Age

Code

# age
mod13 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Age - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod13)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-889.5226  1779.0452  1795.0452  1830.6742  1795.2752   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0085  0.0920     90     no  FocalSpL 
sigma^2.2  0.1498  0.3870    116     no     RecNo 
sigma^2.3  0.6911  0.8313    640     no    Obs_ID 

Test of Moderators (coefficients 1:5):
F(df1 = 5, df2 = 635) = 8.2299, p-val < .0001

Model Results:

      estimate      se    tval   df    pval    ci.lb   ci.ub      
Age     1.1117  0.8092  1.3739  635  0.1699  -0.4772  2.7007      
AgeA    0.4207  0.0677  6.2134  635  <.0001   0.2877  0.5536  *** 
AgeE    0.3229  0.3452  0.9356  635  0.3498  -0.3548  1.0007      
AgeJ    0.2881  0.3484  0.8269  635  0.4086  -0.3961  0.9723      
AgeN    0.3266  0.1911  1.7093  635  0.0879  -0.0486  0.7019    . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod13)*100, 2)

   R2_marginal R2_conditional 
          0.33          18.90 

Code

orchard_plot(mod13,
             mod = "Age",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Predactor type

Code

# type of predator

dat$PredTo[dat$PredTo == ""] <- NA
mod14 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ PredTo - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod14)

Multivariate Meta-Analysis Model (k = 621; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-861.2315  1722.4629  1734.4629  1761.0219  1734.6004   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0003     85     no  FocalSpL 
sigma^2.2  0.1572  0.3965    110     no     RecNo 
sigma^2.3  0.6805  0.8249    621     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 618) = 12.3559, p-val < .0001

Model Results:

         estimate      se    tval   df    pval    ci.lb   ci.ub      
PredToA    0.4284  0.0814  5.2628  618  <.0001   0.2686  0.5883  *** 
PredToB    0.3772  0.3065  1.2306  618  0.2190  -0.2247  0.9791      
PredToN    0.3337  0.0915  3.6466  618  0.0003   0.1540  0.5135  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod14)*100, 2)

   R2_marginal R2_conditional 
          0.25          18.97 

Code

orchard_plot(mod14,
             mod = "PredTo",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Treatment duration

Code

# treatment duration

dat$ln_duration <- log(dat$duration_days)

mod15 <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ln_duration,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod15)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-889.1410  1778.2821  1788.2821  1810.5738  1788.3770   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     90     no  FocalSpL 
sigma^2.2  0.1518  0.3897    116     no     RecNo 
sigma^2.3  0.6888  0.8300    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 8.4146, p-val = 0.0039

Model Results:

             estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt        0.3710  0.0634   5.8526  638  <.0001   0.2465   0.4955  *** 
ln_duration   -0.0455  0.0157  -2.9008  638  0.0039  -0.0763  -0.0147   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod16 <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ln_duration*Type,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod16)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-880.1624  1760.3248  1778.3248  1818.3933  1778.6133   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     90     no  FocalSpL 
sigma^2.2  0.1427  0.3777    116     no     RecNo 
sigma^2.3  0.6882  0.8296    640     no    Obs_ID 

Test of Moderators (coefficients 2:6):
F(df1 = 5, df2 = 634) = 3.5544, p-val = 0.0035

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
intrcpt                        0.4394  0.0717   6.1286  634  <.0001   0.2986 
ln_duration                   -0.0411  0.0173  -2.3802  634  0.0176  -0.0750 
TypeLifeHistory               -0.3341  0.3591  -0.9302  634  0.3526  -1.0393 
TypePhysiology                -0.4638  0.1637  -2.8326  634  0.0048  -0.7853 
ln_duration:TypeLifeHistory    0.0837  0.1012   0.8265  634  0.4088  -0.1151 
ln_duration:TypePhysiology     0.0438  0.0495   0.8849  634  0.3765  -0.0534 
                               ci.ub      
intrcpt                       0.5802  *** 
ln_duration                  -0.0072    * 
TypeLifeHistory               0.3712      
TypePhysiology               -0.1423   ** 
ln_duration:TypeLifeHistory   0.2825      
ln_duration:TypePhysiology    0.1409      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod15)*100, 2)

   R2_marginal R2_conditional 
          3.43          20.87 

Code

bubble_plot(mod15,
             mod = "ln_duration",
             group = "RecNo",
             xlab = "log(duration in days)",
             g = TRUE) +
    geom_point(data = dat,
    aes(x = ln_duration, y = SMD,
    color = Type,
    fill = Type,
    size = 1/sqrt(Vd)), alpha = 0.6) +
    scale_color_discrete() + #+ # how to put the legend for colour
    guides(color = "legend")

Code

#p + geom_point(aes(colour = Type))
#scale_colour_manual(values = c("red", "blue", "green"))

8 Meta-regression: multi-moderator

Code

#######################
# Mulit-variable models
#######################

dat_short$sln_duration <- scale(log(dat_short$duration_days))

mod_full <- rma.mv(yi = SMD, 
               V = VCVs, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ #Design +
                         sln_duration*Type +
                         sln_duration*Treat_mod +
                         Sex,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat_short)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-826.1989  1652.3977  1682.3977  1748.0486  1683.2369   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     81     no  FocalSpL 
sigma^2.2  0.1540  0.3925    104     no     RecNo 
sigma^2.3  0.7169  0.8467    600     no    Obs_ID 

Test of Moderators (coefficients 2:12):
F(df1 = 11, df2 = 588) = 2.0578, p-val = 0.0216

Model Results:

                              estimate      se     tval   df    pval    ci.lb 
intrcpt                         0.4853  0.1076   4.5122  588  <.0001   0.2741 
sln_duration                   -0.1333  0.0896  -1.4886  588  0.1371  -0.3093 
TypeLifeHistory                -0.0267  0.4683  -0.0569  588  0.9546  -0.9463 
TypePhysiology                 -0.5076  0.1840  -2.7580  588  0.0060  -0.8691 
Treat_modV                      0.0477  0.1461   0.3264  588  0.7442  -0.2393 
Treat_modAV                     0.1566  0.1656   0.9455  588  0.3448  -0.1687 
SexF                           -0.2748  0.1295  -2.1219  588  0.0343  -0.5291 
SexM                           -0.0870  0.1531  -0.5679  588  0.5703  -0.3877 
sln_duration:TypeLifeHistory    0.1502  0.4474   0.3358  588  0.7372  -0.7284 
sln_duration:TypePhysiology     0.2185  0.1995   1.0952  588  0.2739  -0.1733 
sln_duration:Treat_modV        -0.0308  0.1418  -0.2168  588  0.8284  -0.3093 
sln_duration:Treat_modAV       -0.1650  0.1804  -0.9150  588  0.3606  -0.5193 
                                ci.ub      
intrcpt                        0.6966  *** 
sln_duration                   0.0426      
TypeLifeHistory                0.8930      
TypePhysiology                -0.1461   ** 
Treat_modV                     0.3347      
Treat_modAV                    0.4818      
SexF                          -0.0205    * 
SexM                           0.2138      
sln_duration:TypeLifeHistory   1.0289      
sln_duration:TypePhysiology    0.6103      
sln_duration:Treat_modV        0.2478      
sln_duration:Treat_modAV       0.1892      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_full)*100, 2)

   R2_marginal R2_conditional 
          7.06          23.50 

Code

orchard_plot(mod_full,
             mod = "Type",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

orchard_plot(mod_full,
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

int_type <- mod_results(mod_full, mod = "sln_duration", group = "RecNo", weights = "prop",
                                   by = "Type")

bubble_plot(int_type, group = "RecNo", mod = "sln_duration", xlab = "ln(duration in days)",
                     legend.pos = "top.left", condition.nrow = 3)

Code

int_trt <- mod_results(mod_full, mod = "sln_duration", group = "RecNo", weights = "prop",
                        by = "Treat_mod")

bubble_plot(int_trt, group = "RecNo", mod = "sln_duration", xlab = "ln(duration in days)",
            legend.pos = "top.left", condition.nrow = 3)

Code

# mulit-model selection
candidates <- dredge(mod_full, trace = 2)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |===================================================================== |  98%

Code

# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 5) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 5)) 

# relative importance of each predictor for all the models
importance <- sw(candidates)

9 Publication bias

Funnel plot: uni-moderator

Funnel plot: multi-moderator

Code

funnel(mod0r, 
       yaxis="seinv",
       type = "rstudent")

Code

funnel(mod_full, 
       yaxis="seinv",
       type = "rstudent")

Egger regression: uni-moderator

Code

# Egger

dat$effectN <- (dat$Ncontrol * dat$NTreat) / (dat$Ncontrol + dat$NTreat)
dat$sqeffectN <- sqrt(dat$effectN)

mod0e <- rma.mv(yi = SMD,
                V = VCV,
                mods = ~ sqeffectN,
                random = list(#~1 | Phylo,
                  ~1 | FocalSpL,
                  ~1 | RecNo,
                  #~1 | SubjectID, # incoprated as VCV
                  ~1 | Obs_ID),
                #R = list(Phylo = cor_tree),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod0e)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-893.0018  1786.0037  1796.0037  1818.2954  1796.0986   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0003     90     no  FocalSpL 
sigma^2.2  0.1483  0.3850    116     no     RecNo 
sigma^2.3  0.6929  0.8324    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 1.0650, p-val = 0.3025

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub      
intrcpt      0.5044  0.1194   4.2234  638  <.0001   0.2699  0.7390  *** 
sqeffectN   -0.0291  0.0282  -1.0320  638  0.3025  -0.0844  0.0262      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(mod0e,
            mod = "sqeffectN",
            group = "RecNo",
            xlab = "Effective N",
            g = TRUE)

Decline effect: uni-moderator

Code

# decline effect
mod0d <- rma.mv(yi = SMD,
                V = VCV,
                mods = ~ Year,
                random = list(#~1 | Phylo,
                  ~1 | FocalSpL,
                  ~1 | RecNo,
                  #~1 | SubjectID, # incoprated as VCV
                  ~1 | Obs_ID),
                #R = list(Phylo = cor_tree),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod0d)

Multivariate Meta-Analysis Model (k = 640; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-892.5966  1785.1932  1795.1932  1817.4849  1795.2881   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0103  0.1016     90     no  FocalSpL 
sigma^2.2  0.1450  0.3808    116     no     RecNo 
sigma^2.3  0.6886  0.8298    640     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 638) = 1.5062, p-val = 0.2202

Model Results:

         estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt   23.4292  18.7537   1.2493  638  0.2120  -13.3973  60.2556    
Year      -0.0114   0.0093  -1.2273  638  0.2202   -0.0297   0.0069    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(mod0d,
            mod = "Year",
            group = "RecNo",
            xlab = "Publication year",
            g = TRUE)

All together

Code

# full model
dat_short$effectN <- (dat_short$Ncontrol * dat_short$NTreat) / (dat_short$Ncontrol + dat_short$NTreat)
dat_short$sqeffectN <- sqrt(dat_short$effectN)

mod_fulle <- rma.mv(yi = SMD, 
                   V = VCVs, 
                   random = list(~1|FocalSpL , 
                                 ~1 | RecNo,
                                 ~1 | Obs_ID),
                   mods =  ~ sqeffectN +
                     Year +
                     sln_duration*Type +
                     sln_duration*Treat_mod +
                     Sex,
                   test = "t",
                   method = "REML",
                   sparse = TRUE,
                   data = dat_short)
summary(mod_fulle)

Multivariate Meta-Analysis Model (k = 600; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-823.4875  1646.9750  1680.9750  1755.3215  1682.0525   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     81     no  FocalSpL 
sigma^2.2  0.1587  0.3983    104     no     RecNo 
sigma^2.3  0.7183  0.8475    600     no    Obs_ID 

Test of Moderators (coefficients 2:14):
F(df1 = 13, df2 = 586) = 1.7639, p-val = 0.0454

Model Results:

                              estimate       se     tval   df    pval     ci.lb 
intrcpt                        14.6850  22.7802   0.6446  586  0.5194  -30.0557 
sqeffectN                      -0.0007   0.0328  -0.0216  586  0.9827   -0.0651 
Year                           -0.0070   0.0113  -0.6221  586  0.5341   -0.0293 
sln_duration                   -0.1252   0.0936  -1.3370  586  0.1817   -0.3090 
TypeLifeHistory                -0.0173   0.4714  -0.0366  586  0.9708   -0.9431 
TypePhysiology                 -0.5322   0.1894  -2.8094  586  0.0051   -0.9043 
Treat_modV                      0.0368   0.1478   0.2489  586  0.8035   -0.2535 
Treat_modAV                     0.1375   0.1690   0.8136  586  0.4162   -0.1944 
SexF                           -0.2740   0.1300  -2.1076  586  0.0355   -0.5293 
SexM                           -0.0953   0.1543  -0.6178  586  0.5369   -0.3983 
sln_duration:TypeLifeHistory    0.1450   0.4510   0.3216  586  0.7479   -0.7407 
sln_duration:TypePhysiology     0.2438   0.2055   1.1863  586  0.2360   -0.1599 
sln_duration:Treat_modV        -0.0261   0.1440  -0.1813  586  0.8562   -0.3088 
sln_duration:Treat_modAV       -0.1497   0.1834  -0.8161  586  0.4148   -0.5099 
                                ci.ub     
intrcpt                       59.4257     
sqeffectN                      0.0637     
Year                           0.0152     
sln_duration                   0.0587     
TypeLifeHistory                0.9086     
TypePhysiology                -0.1602  ** 
Treat_modV                     0.3271     
Treat_modAV                    0.4693     
SexF                          -0.0187   * 
SexM                           0.2077     
sln_duration:TypeLifeHistory   1.0307     
sln_duration:TypePhysiology    0.6475     
sln_duration:Treat_modV        0.2566     
sln_duration:Treat_modAV       0.2105     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

dat_fulle <- qdrg(object = mod_fulle, 
                  data = dat_short)
# marginalized overall mean at vi = 0 and year.c = 0; also weights = "prop" or "cells" average things over proportionally. if not specified, all groups (levels) get the same weights
# res_fulle1 <- emmeans(dat_fulle, 
#                      specs = ~ sqeffectN,
#                      df = mod_fulle$ddf, 
#                      weights = "prop")

10 R Session Informtion

Code

# pander for making it look nicer
sessionInfo() %>% pander()

R version 4.2.3 (2023-03-15)

Platform: aarch64-apple-darwin20 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: pander(v.0.6.5), kableExtra(v.1.3.4), MuMIn(v.1.47.5), emmeans(v.1.8.6), clubSandwich(v.0.5.8), ape(v.5.7-1), easyalluvial(v.0.3.1), ggalluvial(v.0.12.5), alluvial(v.0.1-2), patchwork(v.1.1.2), metafor(v.4.2-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), orchaRd(v.2.0), lme4(v.1.1-32), Matrix(v.1.5-3), here(v.1.0.1), lubridate(v.1.9.2), forcats(v.1.0.0), stringr(v.1.5.0), dplyr(v.1.1.2), purrr(v.1.0.1), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.2) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): TH.data(v.1.1-2), minqa(v.1.2.5), colorspace(v.2.1-0), class(v.7.3-21), ggridges(v.0.5.4), rprojroot(v.2.0.3), estimability(v.1.4.1), rstudioapi(v.0.14), listenv(v.0.9.0), prodlim(v.2023.03.31), fansi(v.1.0.4), mvtnorm(v.1.2-2), mathjaxr(v.1.6-0), xml2(v.1.3.4), codetools(v.0.2-19), splines(v.4.2.3), knitr(v.1.43), jsonlite(v.1.8.5), nloptr(v.2.0.3), compiler(v.4.2.3), httr(v.1.4.6), fastmap(v.1.1.1), cli(v.3.6.1), htmltools(v.0.5.5), tools(v.4.2.3), coda(v.0.19-4), gtable(v.0.3.3), glue(v.1.6.2), Rcpp(v.1.0.10), vctrs(v.0.6.2), svglite(v.2.1.1), nlme(v.3.1-162), progressr(v.0.13.0), timeDate(v.4022.108), gower(v.1.0.1), xfun(v.0.39), globals(v.0.16.2), rvest(v.1.0.3), timechange(v.0.2.0), lifecycle(v.1.0.3), future(v.1.32.0), MASS(v.7.3-58.2), zoo(v.1.8-12), scales(v.1.2.1), ipred(v.0.9-14), hms(v.1.1.3), parallel(v.4.2.3), sandwich(v.3.0-2), RColorBrewer(v.1.1-3), yaml(v.2.3.7), gridExtra(v.2.3), rpart(v.4.1.19), stringi(v.1.7.12), hardhat(v.1.3.0), boot(v.1.3-28.1), lava(v.1.7.2.1), rlang(v.1.1.1), pkgconfig(v.2.0.3), systemfonts(v.1.0.4), evaluate(v.0.21), lattice(v.0.20-45), recipes(v.1.0.6), htmlwidgets(v.1.6.2), tidyselect(v.1.2.0), parallelly(v.1.35.0), magrittr(v.2.0.3), R6(v.2.5.1), generics(v.0.1.3), multcomp(v.1.4-23), pillar(v.1.9.0), withr(v.2.5.0), survival(v.3.5-3), nnet(v.7.3-18), future.apply(v.1.11.0), utf8(v.1.2.3), tzdb(v.0.4.0), rmarkdown(v.2.22), grid(v.4.2.3), data.table(v.1.14.8), webshot(v.0.5.4), digest(v.0.6.31), xtable(v.1.8-4), stats4(v.4.2.3), munsell(v.0.5.0) and viridisLite(v.0.4.2)